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Abstract. Modelling incompressible ideal fluids as a finite collection of 
vortex filaments is important in physics (super-fluidity, models for the 
onset of turbulence) as well as for numerical algorithms used in com- 
puter graphics for the real time simulation of smoke. Here we introduce 
a time-discrete evolution equation for arbitrary closed polygons in 3- 
space that is a discretisation of the localised induction approximation of 
filament motion. This discretisation shares with its continuum limit the 
property that it is a completely integrable system. We apply this poly- 
gon evolution to a significant improvement of the numerical algorithms 
used in Computer Graphics. 



1. Introduction 

The motion of vortex filaments in an incompressible, inviscid fluid has 
aroused considerable interest in quite different areas: 

Differential geometry. The limiting case of infinitely thin vortex filaments 
leads to an evolution equation for closed space curves 7, 

7 = y x 7". (i) 

Equation ((TJ) was discovered in the beginning of the 20th century by Levi- 
Civita and his student Da Rios pQ and is called the smoke ring flow or 
localised induction approximation. In 1972 Hasimoto [2] discovered that 
the smoke ring flow is in fact a completely integrable Hamiltonian system 
equivalent to the non-linear Schrodinger equation. See [3] for more details 
on the history of the smoke ring equation. Subsequently the smoke ring flow 
has been studied by differential geometers as a natural evolution equation 
for space curves [U [6j [7] . Also discrete versions of the smoke ring flow 
in the form of completely integrable evolution equations for polygons with 
fixed edge length have been developed [8j [91 [TO]. 

Fluid dynamics. As will be explained below, for applications in fluid 
mechanics a finite thickness of the vortex filaments has to be taken into 
account. The transition from infinitely thin filaments to filaments of finite 
thickness involves the incorporation of long range interactions (governed by 
the Biot-Savart law) between different filaments and different parts of the 
same filament into the purely local evolution equation ([T]). The resulting 
evolution of vortex filaments has been extensively studied both numerically 
and in the context of explaining the onset of turbulence [llj . Including 
in addition a small amount of viscosity in the equations leads to striking 
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physical effects like vortex reconnection \12\ [T3| [H] and numerical techniques 
like "hairpin removal" \15\ [T6] . 



Computer graphics. Filament-based methods for fluid simulation are 
becoming important in Computer Graphics for special effects in motion pic- 
tures and for real time applications like computer games [171 [T8] . Here the 
emphasis is on physical correctness and speed rather than numerical accu- 
racy. Filament methods are ideal for these applications because complicated 
fluid motions can be created by a graphics designer by modelling the ini- 
tial positions and strengths of the filaments. Moreover, filament methods 
work in unbounded space rather than in a bounded box (as is the case for 
grid-based methods [H]). This is desirable for the simulation of smoke. 



The main goal of this paper is to improve the numerical algorithms cur- 
rently used in Computer Graphics by applying the recent knowledge from 
Discrete Differential Geometry to the motion of polygonal smoke rings. Our 
method makes it possible to model thin filaments by polygons with arbi- 
trarily few vertices. For comparison, using current methods to model a 
circular smoke ring which is thin enough to entrain smoke in a torus shape, 
it necessary to use a regular polygon with at least 800 vertices. 

In Section [2] we will explain the evolution equation for systems of vortex 
filaments that we will discretise. The resulting equations of motion are still 
Hamiltonian like the smoke ring flow (pQ). However, since already Poincare 
knew that a system of vortex filaments consisting of more than three parallel 
straight lines (the "n- vortex problem") fails to be an integrable system [20|. 
p. 58f], we do not believe that this system is an integrable Hamiltonian 
system. Nevertheless it is a small perturbation of the integrable system 
constituted by the limit of infinitely thin filaments. This might be interesting 
for future investigations along the lines of KAM theory. 

In Section [3] we consider polygonal vortex filaments. In this case, there is 
an elementary formula (llip for the Biot-Savart integral. 

In Section0]we will develop an extension of the known discrete-time smoke 
ring flow for polygons of constant edge lengths to arbitrary polygons. This 
is needed because after including the long range Biot-Savart interactions, 
the lengths of the edges will be no longer constant in time. 

In the theory of integrable systems it is known at least since the 1980s 
that integrable difference equations may be interpreted as Darboux trans- 
formations of integrable differential equations [21| , I22 | |2*5]. In the meantime, 
this seminal discovery has lead to a reversed point of view, where the discrete 
integrable systems are considered fundamental and the continuous systems 
appear as smooth limits (see for example |24] and the references therein) . In 
this vein, we will in Section 2] define the discrete-time integrable system in 
terms of iterated Darboux transformations of polygons and show afterwards 
that the smoke ring flow is obtained as a smooth limit. 

In Section [5] we will describe our numerical method that very efficiently 
models the motion of fluids near the smoke ring limit. 
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2. Euler's Equation for Vortex Filaments 

Consider an incompressible, inviscid fluid in euclidean 3-space whose ve- 
locity field u vanishes at infinity and whose vorticity uj = curl u is compactly 
supported. Then u can be reconstructed from uj by the Biot-Savart formula 

f x — z 

u(x) = — — / ^xuj(z)dz. (2) 

47r Jr3 ||a; - z\\ 

The equation of motion can then be written as 

to = [uj, u], (3) 

Viewed as an evolution equation on the vector space M of compactly sup- 
ported divergence-free vector fields on M 3 this is a Hamiltonian system: A 
symplectic form o on M. is defined as follows. Let uj € j\A and u, uj G T^Ad. 
Then 

a L0 {uj,uj) = I det (uj,lj,uj). (4) 
Jk. 3 

Let H : M. — > E be the quadratic function 

H = JJ {U \b-vf dxdy ' (5) 

where (•, •) is the euclidean scalar product on M 3 . Then H is the Hamil- 
tonian for the dynamical system See |2Ut [23] for more details on this 
Hamiltonian description of ideal fluids. 

If the vorticity of a fluid is concentrated on a closed curve 7 in a delta- 
function like manner, by Equation ^) the resulting velocity field u becomes 

" w = -£/F^ xy<s)<is - (6> 

Here T is the circulation around the filament. The problem with Equa- 
tion (|6|) is that in order to determine the motion of 7 itself, u has to be 
evaluated on 7, which results in a logarithmically divergent integral. Usu- 
ally, this problem is handled by considering a vorticity field concentrated 
in a tube around 7 of small but finite radius r. For small r the velocity 
in this tube is dominated by a term proportional to the localised induction 
approximation. (See, for example, [26j p. 36f].) Here we want to derive 
the smoke ring flow by taking the limit r — > 0. In order to prevent vortex 
filaments acquiring infinite speed, one has to scale the circulation T down 
to zero when performing the limit to infinitely thin filaments. This means 
that the fluid velocity ([6]) goes to zero as well. 

The resulting picture is then as follows: The fluid is completely at rest 
away from the filaments while the filaments just cut through the fluid with 
finite speed according to the smoke ring flow: 

l\rjX~'r (7) 
Here the constants Kj account for the fact that the circulation of the differ- 
ent filaments jj might go to zero at a different rate. 

Equation (|7|) can be viewed as a completely integrable Hamiltonian system 
on the space of weighted links (see Figure [T]) endowed with the symplectic 
form 
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^(7,7) = ^^j / det( 7 -,7,7)- (8) 
3 J *H 

For single curves this symplectic form is due to V. I. Arnold [2D]. The 
corresponding Hamiltonian is a weighted sum of the filament lengths 

H = Y,K j Length( 7i ). 

3 

Equation ([7]) can be obtained (using a simple renormalisation of time) as 
a limit as a — > of the following system: Stick with ([8]) as the symplectic 
form, with Kj replaced by the non-zero circulation Tj around jj. As a 
Hamiltonian, use 



hi 87F \J a2 + Il7i(s) -7j(»)ll 

The resulting equation of motion is 

*w -Erf m^m x y (5) (9) 

Y 4 ^ 7 V« 2 + l7,(^)-7 J (5)l 2 
This equation of motion ([DJ can also be derived as follows: 

• Smooth the delta-function like vorticity field ujq of the link by a 
suitable convolution kernel and obtain 

/ \ 3a 2 /" u) Q (y) 
uj[x) = — — / — ,, ay. 

4vr J R3 y / a 2 + \ x _ y \2 b 

• Compute the corresponding velocity field u with curlu = u>: 

<fr) = -£E/-~ =fg^=p »Tt«*- («>) 

* n j J ^Ja l + \x — 7j(s)r 

• Evaluate u on the filaments to obtain ([9]). 
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To summarise: We model fluid motion near the filament limit by a Hamilton- 
ian system on the space of weighted links. This system is still Hamiltonian 
but no longer integrable. Nevertheless it still has all the constants of motion 
induced by invariance with respect to the euclidean symmetry group. For 
example the weighted sum of the area vectors 

3 

is one of the preserved quantities. (Compare Theorem 0] of Section 01) 

The physical approximation implicit in our model is that we ignore pos- 
sible deformations of the internal structure of the filaments and reduce ev- 
erything to the evolution of the filament curves. The finite thickness of the 
filaments is taken into account by applying a fixed convolution kernel. 



3. Polygonal Vortex Filaments 

In order to develop a numerical method for modelling fluid motion near 
the filament limit we have to discretise the vortex filaments, i.e. we replace 
them by polygons. If 7 is a piecewise linear parametrisation of a closed 
polygon, on each edge we have 7" = and we find an explicit anti-derivative 
for the integrands of equation (fT0|) : 

(7 ' Y) 7 X7'V= 3 - (11) 

y^ + M 2 (l7'l 2 a2 + | 7 X7'|2) ') y^TW 

Here we have abbreviated x — Jj(s) to 7, 7j(s) to 7' and the prime is deriva- 
tion with respect to s. 

Inspection of Equation (j lip reveals the following problem: The two adja- 
cent edges have no influence at all on the velocity of a vertex. This amounts 
to effectively employing a distance cut-off in order to regularise the singular 
integral ([6]) for points on 7. It is known [26] that this is roughly equivalent to 
modelling vortex tubes of thickness equal to the edge length of the polygon. 
Using this model we would therefore be unable to model thin (and therefore 
fast) filaments without using excessively many edges for each polygon. 

The contribution of local effects behaves like the smoke ring flow and the 
resulting equation of motion for a vertex 7 j of a polygonal vortex filament 
7 is then 

7i = u(ji) + \Kibi, (12) 

where u is given by Equation (|1U|) using (|11|) . Kj6j denotes curvature times 
binomial at ji, and A is constant for fixed a. Since the non-local effects 
quickly destroy any arc-length parametrisation (i.e. the lengths of the dif- 
ferent edges of the polygon) and we do not have an adequate notion of 
curvature for arbitrary polygons, we can not evaluate (fl2|) directly. 

On the other hand, for polygons with constant edge lengths it is known 
that the doubly discrete smoke ring (or Hasimoto) flow [9] captures excel- 
lently the qualitative behaviour of the smooth smoke ring flow. In the next 
section we will discuss a version of this doubly discrete smoke ring flow which 
works also for polygons with varying edge lengths. 
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Figure 2. A polygon 7 and an edge of its Darboux trans- 
form 77. 

4. Darboux Transformation of Polygons 

In this section we develop a discrete-time evolution for closed polygons 
that has the smoke ring flow (fTJ) as a limit when the polygon approaches 
a smooth curve and the time-step goes to zero. This evolution (obtained 
by iterating so-called Darboux transformations) shares with its continuum 
limit the property that it is a completely integrable system in the sense 
that it comes from a Lax pair of quaternionic 2 x 2-matrices with a spectral 
parameter. (This system therefore fits into the framework of [27].) The 
constants of the motion of the discrete system converge to constants of the 
motion of the smooth system in the limit. 

Let 7 : Z — > ]R 3 be an immersed polygon in M 3 , where immersed means 
that 7i 7^ 7i+i for all j£Z, and let Si = ji+i — Ji- If 7 is periodic with some 
period n, then the polygon is closed and 7 may be interpreted as a function 
on Z/nZ. In the following, we identify IR 3 with the imaginary quaternions 
ImH = {xi + yj + zk \ x, y, z G M}. 

Definition. A polygon r/ is called a Darboux transform of 7 with twist 
parameter r 6 ]R and distance / > 0, if 1 1 77^ — 7i|| = I for all i £ Z, and 
the normalised difference vectors Tj defined by ITi = rji — 7j satisfy the 
quaternionic equation 



The Darboux transformation of polygons and its relationship with the 
nonlinear Schrodinger equation and smoke ring flow was treated in [9] under 
the assumption that the polygon 7 has constant edge length. To drop this 
assumption was suggested to us by Tim Hoffmann |28j . 

Geometrically, Equation (|13p has the following meaning (see Figure [2|) . 
The difference vector Tj+i is obtained from Tj by a rotation with axis ITi — 
Si. The quadrilateral 7i7i+i 7 ?i+i 7 ?« is therefore a "folded parallelogram". In 
particular, corresponding edges of 7 and 77 have the same length. The angle 
of rotation is 2 arctan(||/Tj — Si\\/r). For r = it is it. For r —> ±00, it goes 
to zero and in the limit the Darboux transformation becomes a translation. 

Equation f) 131) can be written in the form 



T i+ i = (-r + m - Si)Ti{-r + lT t - Si) 



(13) 



T i+1 = (aTi + b)(cTi + d) 



-1 



(14) 
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where a,b,c,d G EI depend on Si and the parameters l,r. That is, for 
each i £ Z, T l+ \ is obtained by applying a quaternionic fractional linear 
transformation /, : I ^ H to T u where i = HU {oo}. Indeed, (USD is 
equivalent to 

T i+1 = (ITi -r-St) ((r + S^I- + Z) _1 . (15) 

To see this note that = —Ti because T, is a purely imaginary unit 
quaternion, and hence Tj(— r + ZT, — S^) -1 = (rT, + I + SjTj). 

It is convenient to rewrite fractional linear transformations as matrix mul- 
tiplication. Just as the extended complex plane C = CU {00} can be identi- 
fied with the Riemann sphere S 2 and with the complex projective line CP 1 , 
H = S 4 = HP 1 . The quaternionic projective line HP 1 is the set of (quater- 
nionic) 1-dimensional subspaces of the vector space H 2 over HI. We consider 
H 2 as right vector space: the product of a vector ( p ) G H 2 and a scalar 

ASM is K)A = (g). A point 



G MP 1 




corresponds to the point pq~ l G H, and p, q are quaternionic homogeneous 
coordinates for this point. Now any fractional linear transformation of H can 
be written as quaternionic 2 x 2-matrix acting from the left on quaternionic 
homogeneous coordinates of HP 1 : Writing Ti in quaternionic homogeneous 
coordinates, 

Ti = T^(TP)-\ 

one obtains from f) 1 51) 

(r|) =Wr) (^)' m,p):= (p + ^ " P A Si )- (16) 

The following Theorem [T] characterises the Darboux transformations of 
polygons via a Lax pair of quaternionic 2 x 2-matrices with spectral param- 
eter. Theorem [2] is a permutability theorem for these Darboux transforma- 
tions. 

Theorem 1 (Lax pair). Let Si = 7^+1 — ji, |Tj| = 1, and let Ui(X, p) be 
defined by 1116}) and 



Vi(X,p) = 



X -p-Si 
p + Si x 



X -p + r - ITi 

r + m X 



Then 



V i+1 (X,p)Ui(X,p) = Ui(X,p)Vi(X,p) (17) 
for all A, p G R, if and only if S and T satisfy (T3\) and 

lT l+1 + Si = Si + lTi. (18) 

That is, if and only if r] = 7 + IT is a Darboux transform of 7 with twist 
parameter r and distance I, and Si = rji + i — rfc. 
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Of course (|1T|) means that the following diagram commutes: 

"v i+1 

B 2 

Proof. Note that in general for quaternionic 2 x 2-matrices with a+, a, b, b £ 
H and A € M the equality 

/ A a + A / A b A ( A 6 A / A a A 
V -a+ X ) \ -b X ) ~ \ -b X ) \ -a X ) 

is equivalent to 

a + 6 = ba and A(a+ + b) = X(b + a). 
It follows that flUD holds for all A € R, if and only if (JTHJ) holds and 

(-p + r - lT i+1 )(-p - Si) = (-p - Si)(-p + r - ITJ. 

Use (I18p to eliminate Si from this equation and gather terms of equal power 
in p on both sides. The coefficients of p 2 are both 1, and the coefficients of 
p are obviously equal. What remains is the equation 

(r - lT i+l ){-Si) = (-Si - IT 1+1 + m)(r - Ity. 

Solve for Tj + i to obtain (|13p . □ 

Theorem 2 (Permutability). Suppose rj = 7 + IT is a Darboux transform 
of 7 with twist parameter r and distance I, and fj = 7 + AT is a Darboux 
transform of 7 with twist parameter p and distance A, then r/ + AT with 

f= (XT - p + r - lT)({p - r + lT)f + X)' 1 (19) 

is a Darboux transformation of rj with twist parameter p and distance A. 

Proof. Note that T is obtained by applying the quaternionic fractional linear 
transformation represented by the matrix Vi(X, p) to Tj. Let us write T = 
Vi(X, p)Ti for short. Since r) is a Darboux transform of 7 with twist parameter 
p and distance A, Equation (fT6|) says that T + i = C/j(A, /3)Tj. Now TheoremQ] 
implies T + i = Ui(X, p)% and hence (again by Equation ([To]) ). 7/ + AT is a 
Darboux transformation of rj with twist parameter p and distance A. □ 

Even if 7 is a closed curve, the curves obtained by iterating (| 13j) will 
in general not close up. However, we will see that any closed curve has 
generically two closed Darboux transforms. 

The fractional linear transformations : Tj 1— > T l+ i that are represented 
by the matrices Ui(l,r) have the special property that they map the unit 
sphere S 2 = {q G ImH | q 2 = —1} to itself. This follows directly from 
Hence the restrictions fi\$2 are Mobius transformations of S 2 . In fact, 
they are orientation preserving Mobius transformations: By continuity, it is 
enough to check this for a particular value of r and I; and for r = 0, I = 
one obtains T+i = SiTiS^ 1 , which is a 180° rotation with axis Si. 

In order to find for given /, r the closed Darboux transforms of 7, one 
has to look for choices of the initial unit vector Tq such that the recursion 
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(fT5|) generates a sequence with period n, i.e. To = T n . The composition 
fn-i ° • • • ° fo, which maps To l— * T n , is represented by the monodromy 
matrix 

Hi, r = U n - 1 (l,r)---U 2 (l,r)Ui(l,r)U {l,r). 

It is is itself an orientation-preserving Mobius transformation of the unit 
sphere S 2 onto itself. For special cases (we will see below that this cannot 
happen for all I, r) this Mobius-transformation could be the identity, but in 
general it will have exactly two fixed points (counted with multiplicity). 
With each closed curve 7 we have thus associated a monodromy map 

fn-i ■ ■ ■ fo- To will be a fixed point if and only if (~i J is an eigenvector 
of the monodromy matrix H^ r . The following theorem is an immediate 
consequence of Theorem [TJ 

Theorem 3. Suppose r] = 7 + IT is a closed Darboux transform of 7 with 
distance I and twist parameter r. Then for all X and p, the monodromy 
matrix of i] is conjugate to the monodromy matrix H\ p 0/7: 

Hl p = V (X,p)H XtP V (X,p)- 1 . (20) 

This means that if ("1°) *s an eigenvector of H\ p , then Vo(A,p) ^ s an 

eigenvector of H\ p - 

Moreover, one can compute all closed Darboux transforms of rj without 
having to solve an eigenvalue problem, even without iterating the /j. Indeed, 
by Theorem [21 all closed Darboux transforms of rj are given by (I19p . 

Theorem [3] implies that apart from the edge lengths there are many other 
quantities connected with closed polygons that are invariant under Darboux 
transforms: For each A, p the conjugacy class of the monodromy matrix 
H\ tP is invariant. We will show that this implies a nice geometric invariant: 
The area vector of a closed polygon turns out to be invariant under Darboux 
transformations (Theorem H]). 

To derive the invariance of the area vector from the invariance of the 
conjugacy class of the monodromy matrix, we equip the set of quaternionic 
2 x 2-matrices of the form 

with the structure of a C-algebra that is isomorphic to gl(2,C). First note 
that a quaternionic 2 x 2-matrix is of the form (I2ip precisely if it commutes 
with 

' -1 
1 

Define the multiplication of such a matrix with a scalar A + ip 6 C by 

(A + ip)A = (XI + pJ)A, (22) 

where / is the identity matrix. 

The complex multiples of the identity are then 

Z = (X + ip)I = XI + pJ= ( X p (23) 



J 
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Thus we can write Ui(X, p) and Vi(X,p) as 

' Si 



Ui(X,p) = (X + ip)I + J 



Si 



Vi(X,p) = {X + ip)I + J 



-r + ITi 
-r + ITi 



Remark. This means we can combine A and p into one complex spectral 
parameter A + ip. 

Equation (|23p also implies that the trace-free complex matrices in gl{2, C) 
correspond to those matrices of the form (|21f) with a,b £ ImH. Further, a 
matrix of the form pip has a, b € ImH precisely if its square is a matrix 
of the form (|23|) . that is, a (complex) multiple (with multiplication defied 
by (I22p ) of the identity. Identifying C with the matrices of the form (|23p we 
obtain 

-tr c ( ? _b ] = Re a + (Re 6) J 



2 \ b a 
and 

detc ( b ~a ) =\({^ A f -^ A2 ) = \ a2 \-\ h2 \ + 2 ^ b ) J - 



In particular 



I -r - S < 



det c(^ r + (S z ) = r - r- 4 - |5r + 2irJ, 

which vanishes precisely when r = 0, / = ±|5|. Using the notation 

S 



diag(S) :- 



5 



for 5* G EI we can express H\ >p as 

= (Z + J diag(S n -i)) ■ ■ ■ (Z + Jdiag(S )), 

with Z given by (I23p . Hence detc is a complex polynomial of degree 
2n with zeroes precisely at Z = ±|5o|, i^n-ij. By Theorem [3] this 
determinant is invariant under Darboux transforms. This just corresponds 
to the fact that the edge lengths are invariant by construction. Non-trivial 
further invariants come from the complex polynomial 

P(Z) = tr c H z 

of degree n. Let us look at the polynomial coefficients of Hz itself: 



H z = zkA 



n—k i 

k=l 



where 



In particular, 



A k = J k ^ diag(S jl ■ ■ ■ S jk ] 
n-l>ji,...,j k >0 



A =1 

rk 



M =J k E n k= Mag(S k ) = 0, 
M =~ En-i>i>j>o diag(SiSj). 
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That is, A2 is a diagonal matrix with both diagonal entries equal to 

q = - S i S r 
n—l>i>j>0 

The real part of q is 

R e (<?) = En-l>i>j>o(^' ^j) = \ Y^ijtji^ij &j) = \\ Ei=0 Si\ 2 — \ Ei=o \Si\ 

_ 1 sr^n-l I q 12 

This is a function of the edge lengths and therefore not interesting. The 
imaginary part of q is given by 

2A := lm(q) = - £ 4> . Si x Sj 

= E"=i (£?=i Si) x S, 

= Ep^Tj - To) x ( 7i+ i - 7j) 

= E"=i(7j - 7o) x (7 i+ i - 70). 

This invariant A is just the area vector. The following proposition (with 
obvious proof) clarifies its geometrical meaning. 

Proposition 1. Let a € M 3 &e a unit vector, \a\ = 1, and endow the plane 
a ± with the volume form 

det a x(X, Y) := det M3 (a, X, Y). 

Then the area enclosed by the orthogonal projection 7 of the polygon 7 

In = In ~ (ln,a)a 

is equal to (M, a) . 

This explains the name area vector: It encodes all the projected areas. 
Theorem 4. The area vector A is invariant under Darboux transforms. 



Proof. By (J20J), the monodromy matrix of the Darboux transformed curve 

n 
k=0 

satisfies 

H\{Z + J{-rI + I diagT )) = (Z + J(-rI + I diagT ))H z (24) 

Using 

H z = Z n + Z n - 2 A 2 + ... + A , 
H\ = Z n + Z n ~ 2 A v 2 + ... + A v 

and comparing the Z n-2 -coefficients in both sides of (|24|) we obtain A\ = A 2 . 

□ 

Finally we consider the continuum limit of smooth curves 7 : S 1 — ► M 3 
and indicate why Darboux transforms with small parameters I, r do indeed 
converge to the smoke ring flow ([!]). The continuum limit of (|13|) is obtained 
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by replacing S by hS and then computing T' := -£i\ h=0 Th- The resulting 
differential equation is 

T' = (TS - ST^-r + lTy 1 

or 

T' = x(lTxS- rS), (25) 

r z _|_ ^ 

where S : M — > M 3 is given by 

y = s. 

One can check that, as expected, the transformed curve rj = 7 + IT satisfies 

\v'\ = 

The monodromy of the ODE (|25p is a Mobius transformation of S 2 that 
generically has exactly two fixed points. Thus, for generic parameters I and 
r a space curve 7 has exactly two closed Darboux transforms. 

Assume now that we have for r = — I a family of such closed Darboux 
transforms rji that depend analytically on I. Then we reparametrise r\i as 

71(a) := m (s -l)= 7 (s - I) + lT{s - I). (26) 

Then 70 = 7 and comparing coefficients of I in the power series expansion 
of (EBl) we obtain 



, = °> ^72 



2 

11 = i x 7" 



a 2 

U / w -.// 



9/ «=o '" c/^ 1=0 

Hence 

72-7o = *V x 7 " + 0(/ 3 ). 
A small time-step At of the smoke ring flow is therefore approximated by a 
Darboux transform with length I given by I 2 = At. 

Remark. In order to eliminate the reparametrising effect of the Darboux 
transforms it is convenient to apply first a Darboux transform with param- 
eters I and — r followed by a reverse Darboux transform with parameters I 
and r. This will cancel out the (first order in t) tangential shift and leave 
only the (second order in t) smoke ring evolution (see [8]). 

5. An ALGORITHM FOR THE REAL TIME SIMULATION OF FLUID FLOW 

Based on the theoretic foundations covered in the previous sections, we 
have implemented the following algorithm for the simulation of fluid flow. 
Our aim was to develop an algorithm which is fast enough to generate re- 
alistic looking computer animations of fluid motion in real time. Figure 
shows a sample screen shot from a simulation which runs smoothly on stan- 
dard hardware. We assume the vorticity is concentrated along a few vortex 
rings, which we represent by closed polygons. Their motion is governed by 
a mixture of the velocity field induced by the polygonal vortex rings via the 
smoothed Biot-Savart formula (jlOp of Section [3l and Darboux transforma- 
tions which approximate a time step of the polygonal smoke ring flow as 
explained in Section |H The rationale behind this scheme is that the veloc- 
ity field induced by an edge of a polygonal vortex filament is zero on that 
edge itself. Thus, the adjacent edges do not contribute to the velocity of a 
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FIGURE 3. 256 2 fluid particles evolving under the influence 
of three polygonal vortex filaments. 

vertex. The Darboux transforms make up for this lack of local interaction. 
The following is a summary description of the algorithm. Details (in par- 
ticular how we set the parameters rj and li of the Darboux transformation) 
are given below. 

input: 

• positions 7^ of the jih. vertex of the ith. polygonal vortex fila- 
ment 7i, where i = 1 . . . m, j = 1 . . . raj. 

• strengths Tj and smoothing (thickness) parameters aj of the 
vortex filaments. 

• positions pi € M 3 of advected particles, where i = 1 . . . k. 

• time-step At . 
loop: 

1 Compute a double Darboux transform rji with parameters =Frj, Zj 
of each polygon 7;. 7^ *- t]iy 

2 Solve jij = ui'jij) for time-step At, where u{x) is the velocity 
field obtained by the smoothed Biot-Savart formula (|10l) . 

3 Update the particle positions pi by solving p. L = u(pi) for time- 
step At. 

In Step 1, we determine the parameters li and rj as follows. The amount of 
smoke ring flow needed to make up for the lack of local interaction depends 
on the thickness en, the number of edges rij and the total length Lj of 
7i. Since we do not know the correct speed for an arbitrary polygon, we 
determine the parameters for the test case of a regular nj-gon with same 
strength, thickness and total length. We choose the parameters in such 
a way that for the regular rij-gon the sum of self-induced velocity from 
the Biot-Savart formula ()10p plus the effect of a double Darboux transform 
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coincides with the analytically known speed U for a circle with same length 

^ = ^fln^-l), (27) 
ZLi \ nai J 

compare [26!, p. 212]. We compute the self- induced speed Ui of the n^-gon 
by evaluating the smoothed Biot-Savart formula ()10p at one vertex for all 
edges of the nj-gon. This speed is slower than Ui because the adjacent edges 
have no influence on a vertex, see Section [3l Now we choose and li such 
that a double Darboux transformation translates the regular nj-gon by a 
distance of (Ui — U) At. A single Darboux transform of the regular nj-gon 
is a translation in binormal direction plus a non-zero rotation about the 
centre axis. The rotation cancels out for a double Darboux transform and is 
therefore arbitrary. We choose the rotation angle to be 2n/rii, which leads 
to the following formulas for li and r^. 

k = \j (Li/rii) 2 + of , n = Oi cot(7r/nj) , 

where we have abbreviated \(JJi — U) At by Oi. 

In Step 2, we use the fourth order Runge-Kutta scheme (RK4) to solve the 
ordinary the differential equation x = u(x) for the time-step At. To advect 
the large number of particles in Step 3 we use second order Runge-Kutta 
(RK2), where we use the two polygon positions after Step 1 and Step 2 as 
intermediate values. To improve performance further, this step is computed 
on the computer's graphics chip (GPGPU). 

Evaluating u(x) via Equation (fTUl) is unproblematic, because the integral 
on the right hand side can be solved explicitly for straight line segments; see 
Equation (fTT|) in Section [3l 
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